DJIA,GSPC,NDX,GDAXI,FCHI,SSEC,SENSEX Analysis

https://github.com/cran/VineCopula - Copula family

setwd('~/Documents/VaR')

prices <- read.csv('./data/StockIndexData.csv')

library(VineCopula)
library(copula)
## 
## Attaching package: 'copula'
## The following object is masked from 'package:VineCopula':
## 
##     fitCopula
library(tseries)

DJIA,GSPC,NDX,GDAXI,FCHI,SSEC,SENSEX

prices.DJIA <- prices[,c('DJIA')]

summary(prices.DJIA)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##    6547   10220   11020   12020   13270   18640
n <- length(prices.DJIA)
rtn.DJIA <- rep(0,n-1)

for(i in seq(n-1)){
  rtn.DJIA[i] <- log(prices.DJIA[i+1]/prices.DJIA[i]) * 100.
}

prices.GSPC <- prices[,c('GSPC')]

summary(prices.GSPC)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   676.5  1133.0  1296.0  1364.0  1478.0  2190.0
n <- length(prices.GSPC)
rtn.GSPC <- rep(0,n-1)

for(i in seq(n-1)){
  rtn.GSPC[i] <- log(prices.GSPC[i+1]/prices.GSPC[i]) * 100.
}

prices.NDX <- prices[,c('NDX')]

summary(prices.NDX)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    NA's 
##   938.5  1569.0  1944.0  2394.0  2954.0  4910.0    1063
n <- length(prices.NDX)
rtn.NDX <- rep(0,n-1)

for(i in seq(n-1)){
  rtn.NDX[i] <- log(prices.NDX[i+1]/prices.NDX[i]) * 100.
}

prices.GDAXI <- prices[,c('GDAXI')]

summary(prices.GDAXI)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    NA's 
##    2203    4934    6172    6461    7583   12370      60
n <- length(prices.GDAXI)
rtn.GDAXI <- rep(0,n-1)

for(i in seq(n-1)){
  rtn.GDAXI[i] <- log(prices.GDAXI[i+1]/prices.GDAXI[i]) * 100.
}

prices.FCHI <- prices[,c('FCHI')]

summary(prices.FCHI)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    NA's 
##    2403    3649    4220    4291    4836    6857      50
n <- length(prices.FCHI)
rtn.FCHI <- rep(0,n-1)

for(i in seq(n-1)){
  rtn.FCHI[i] <- log(prices.FCHI[i+1]/prices.FCHI[i]) * 100.
}

prices.SSEC <- prices[,c('SSEC')]

summary(prices.SSEC)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    NA's 
##    1011    1572    2104    2288    2850    6092     317
n <- length(prices.SSEC)
rtn.SSEC <- rep(0,n-1)

for(i in seq(n-1)){
  rtn.SSEC[i] <- log(prices.SSEC[i+1]/prices.SSEC[i]) * 100.
}

prices.SENSEX <- prices[,c('SENSEX')]

summary(prices.SENSEX)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    NA's 
##    2600    4757   13670   13040   18840   29680     202
n <- length(prices.SENSEX)
rtn.SENSEX <- rep(0,n-1)

for(i in seq(n-1)){
  rtn.SENSEX[i] <- log(prices.SENSEX[i+1]/prices.SENSEX[i]) * 100.
}
rtn.DJIA.ar <- ar(rtn.DJIA, order.max = 1, method = 'ols')
rtn.DJIA.ar.resid.garch <- garch(rtn.DJIA.ar$resid[2:length(rtn.DJIA.ar$resid)],trace=FALSE)

rtn.GSPC.ar <- ar(rtn.GSPC, order.max = 1, method = 'ols')
rtn.GSPC.ar.resid.garch <- garch(rtn.GSPC.ar$resid[2:length(rtn.GSPC.ar$resid)],trace=FALSE)

rtn.NDX.1<-na.omit(rtn.NDX)
rtn.NDX.ar <- ar(rtn.NDX.1, order.max = 1, method = 'ols')
rtn.NDX.ar.resid.garch <- garch(rtn.NDX.ar$resid[2:length(rtn.NDX.ar$resid)],trace=FALSE)

rtn.GDAXI.1<-na.omit(rtn.GDAXI)
rtn.GDAXI.ar <- ar(rtn.GDAXI.1, order.max = 1, method = 'ols')
rtn.GDAXI.ar.resid.garch <- garch(rtn.GDAXI.ar$resid[2:length(rtn.GDAXI.ar$resid)],trace=FALSE)

rtn.FCHI.1<-na.omit(rtn.FCHI)
rtn.FCHI.ar <- ar(rtn.FCHI.1, order.max = 1, method = 'ols')
rtn.FCHI.ar.resid.garch <- garch(rtn.FCHI.ar$resid[2:length(rtn.FCHI.ar$resid)],trace=FALSE)

rtn.SSEC.1<-na.omit(rtn.SSEC)
rtn.SSEC.ar <- ar(rtn.SSEC.1, order.max = 1, method = 'ols')
rtn.SSEC.ar.resid.garch <- garch(rtn.SSEC.ar$resid[2:length(rtn.SSEC.ar$resid)],trace=FALSE)

rtn.SENSEX.1<-na.omit(rtn.SENSEX)
rtn.SENSEX.ar <- ar(rtn.SENSEX.1, order.max = 1, method = 'ols')
rtn.SENSEX.ar.resid.garch <- garch(rtn.SENSEX.ar$resid[2:length(rtn.SENSEX.ar$resid)],trace=FALSE)

rm(prices.DJIA, prices.GSPC, prices.NDX, prices.GDAXI, prices.FCHI,prices.SSEC, prices.SENSEX)
rm(rtn.DJIA, rtn.GSPC, rtn.NDX, rtn.GDAXI, rtn.FCHI,rtn.SSEC, rtn.SENSEX)
rm(rtn.NDX.1, rtn.GDAXI.1, rtn.FCHI.1,rtn.SSEC.1, rtn.SENSEX.1)

Standardize the residual from AR(1)XGARCH(1,1)

std.residual.DJIA <- rtn.DJIA.ar$resid/sqrt(rtn.DJIA.ar.resid.garch$fitted.values)
## Warning in sqrt(rtn.DJIA.ar.resid.garch$fitted.values): NaNs produced
## Warning in rtn.DJIA.ar$resid/sqrt(rtn.DJIA.ar.resid.garch$fitted.values):
## longer object length is not a multiple of shorter object length
std.residual.GSPC <- rtn.GSPC.ar$resid/sqrt(rtn.GSPC.ar.resid.garch$fitted.values)
## Warning in sqrt(rtn.GSPC.ar.resid.garch$fitted.values): NaNs produced
## Warning in rtn.GSPC.ar$resid/sqrt(rtn.GSPC.ar.resid.garch$fitted.values):
## longer object length is not a multiple of shorter object length
std.residual.NDX <- rtn.NDX.ar$resid/sqrt(rtn.NDX.ar.resid.garch$fitted.values)
## Warning in sqrt(rtn.NDX.ar.resid.garch$fitted.values): NaNs produced
## Warning in rtn.NDX.ar$resid/sqrt(rtn.NDX.ar.resid.garch$fitted.values):
## longer object length is not a multiple of shorter object length
std.residual.GDAXI <- rtn.GDAXI.ar$resid/sqrt(rtn.GDAXI.ar.resid.garch$fitted.values)
## Warning in sqrt(rtn.GDAXI.ar.resid.garch$fitted.values): NaNs produced
## Warning in rtn.GDAXI.ar$resid/sqrt(rtn.GDAXI.ar.resid.garch$fitted.values):
## longer object length is not a multiple of shorter object length
std.residual.FCHI <- rtn.FCHI.ar$resid/sqrt(rtn.FCHI.ar.resid.garch$fitted.values)
## Warning in sqrt(rtn.FCHI.ar.resid.garch$fitted.values): NaNs produced
## Warning in rtn.FCHI.ar$resid/sqrt(rtn.FCHI.ar.resid.garch$fitted.values):
## longer object length is not a multiple of shorter object length
std.residual.SSEC <- rtn.SSEC.ar$resid/sqrt(rtn.SSEC.ar.resid.garch$fitted.values)
## Warning in sqrt(rtn.SSEC.ar.resid.garch$fitted.values): NaNs produced
## Warning in rtn.SSEC.ar$resid/sqrt(rtn.SSEC.ar.resid.garch$fitted.values):
## longer object length is not a multiple of shorter object length
std.residual.SENSEX <- rtn.SENSEX.ar$resid/sqrt(rtn.SENSEX.ar.resid.garch$fitted.values)
## Warning in sqrt(rtn.SENSEX.ar.resid.garch$fitted.values): NaNs produced
## Warning in rtn.SENSEX.ar$resid/sqrt(rtn.SENSEX.ar.resid.garch
## $fitted.values): longer object length is not a multiple of shorter object
## length
rtns <- cbind(std.residual.DJIA[1067:4532],
              std.residual.GSPC[1067:4532],
              std.residual.NDX[2:3467],
              std.residual.GDAXI[953:4418],
              std.residual.FCHI[967:4432],
              std.residual.SSEC[658:4123],
              std.residual.SENSEX[674:4139]
              )


rm(std.residual.DJIA,std.residual.GSPC,std.residual.NDX,std.residual.GDAXI ,std.residual.FCHI,
   std.residual.SSEC,std.residual.SENSEX)

cor(rtns,method='pearson')
##             [,1]        [,2]          [,3]        [,4]          [,5]
## [1,]  1.00000000  0.97482107  0.4626361048 -0.02897670 -0.0448450319
## [2,]  0.97482107  1.00000000  0.4877918823 -0.02413000 -0.0454992144
## [3,]  0.46263610  0.48779188  1.0000000000 -0.05170464  0.0005284583
## [4,] -0.02897670 -0.02413000 -0.0517046389  1.00000000 -0.0142933609
## [5,] -0.04484503 -0.04549921  0.0005284583 -0.01429336  1.0000000000
## [6,] -0.02312062 -0.02512210  0.0025580560  0.01929367  0.0158593348
## [7,] -0.01925440 -0.02025158 -0.0048500626  0.01343810 -0.0123316220
##              [,6]         [,7]
## [1,] -0.023120619 -0.019254398
## [2,] -0.025122102 -0.020251582
## [3,]  0.002558056 -0.004850063
## [4,]  0.019293672  0.013438099
## [5,]  0.015859335 -0.012331622
## [6,]  1.000000000  0.013798025
## [7,]  0.013798025  1.000000000
cor(rtns,method='kendall')
##              [,1]         [,2]         [,3]         [,4]         [,5]
## [1,]  1.000000000  0.840936111  0.339267541 -0.018745030 -0.020497615
## [2,]  0.840936111  1.000000000  0.370060343 -0.012876103 -0.020819688
## [3,]  0.339267541  0.370060343  1.000000000 -0.026871801 -0.001039327
## [4,] -0.018745030 -0.012876103 -0.026871801  1.000000000 -0.012353858
## [5,] -0.020497615 -0.020819688 -0.001039327 -0.012353858  1.000000000
## [6,] -0.003618578 -0.004117842 -0.000861471  0.010674214 -0.002154094
## [7,] -0.004118841 -0.003936654 -0.004190783 -0.001609867 -0.020418012
##              [,6]         [,7]
## [1,] -0.003618578 -0.004118841
## [2,] -0.004117842 -0.003936654
## [3,] -0.000861471 -0.004190783
## [4,]  0.010674214 -0.001609867
## [5,] -0.002154094 -0.020418012
## [6,]  1.000000000  0.009381258
## [7,]  0.009381258  1.000000000
cor(rtns,method='spearman')
##              [,1]         [,2]         [,3]         [,4]         [,5]
## [1,]  1.000000000  0.962107436  0.454675399 -0.027960277 -0.030258828
## [2,]  0.962107436  1.000000000  0.484046514 -0.019486288 -0.030636594
## [3,]  0.454675399  0.484046514  1.000000000 -0.040064311 -0.001568855
## [4,] -0.027960277 -0.019486288 -0.040064311  1.000000000 -0.018482587
## [5,] -0.030258828 -0.030636594 -0.001568855 -0.018482587  1.000000000
## [6,] -0.005895352 -0.006491070 -0.001236738  0.015962569 -0.002794678
## [7,] -0.006206004 -0.006164209 -0.006189422 -0.002249416 -0.030567593
##              [,6]         [,7]
## [1,] -0.005895352 -0.006206004
## [2,] -0.006491070 -0.006164209
## [3,] -0.001236738 -0.006189422
## [4,]  0.015962569 -0.002249416
## [5,] -0.002794678 -0.030567593
## [6,]  1.000000000  0.013779081
## [7,]  0.013779081  1.000000000
splom2(rtns)

g <- rtns[,1]
h<-hist(g, breaks=10, density=10, col="lightgray", xlab="Accuracy", main="Overall") 
    xfit<-seq(min(g),max(g),length=40) 
    yfit<-dnorm(xfit,mean=mean(g),sd=sd(g)) 
    yfit <- yfit*diff(h$mids[1:2])*length(g) 
    lines(xfit, yfit, col="black", lwd=2)

h
## $breaks
##  [1] -5 -4 -3 -2 -1  0  1  2  3  4  5
## 
## $counts
##  [1]    1   10   70  366 1183 1401  377   53    3    2
## 
## $density
##  [1] 0.0002885170 0.0028851702 0.0201961916 0.1055972302 0.3413156376
##  [6] 0.4042123485 0.1087709175 0.0152914022 0.0008655511 0.0005770340
## 
## $mids
##  [1] -4.5 -3.5 -2.5 -1.5 -0.5  0.5  1.5  2.5  3.5  4.5
## 
## $xname
## [1] "g"
## 
## $equidist
## [1] TRUE
## 
## attr(,"class")
## [1] "histogram"
g <- rtns[,2]
h<-hist(g, breaks=10, density=10, col="lightgray", xlab="Accuracy", main="Overall") 
    xfit<-seq(min(g),max(g),length=40) 
    yfit<-dnorm(xfit,mean=mean(g),sd=sd(g)) 
    yfit <- yfit*diff(h$mids[1:2])*length(g) 
    lines(xfit, yfit, col="black", lwd=2)

h
## $breaks
##  [1] -5 -4 -3 -2 -1  0  1  2  3  4  5
## 
## $counts
##  [1]    4   14   82  362 1139 1409  398   52    4    2
## 
## $density
##  [1] 0.001154068 0.004039238 0.023658396 0.104443162 0.328620889
##  [6] 0.406520485 0.114829775 0.015002885 0.001154068 0.000577034
## 
## $mids
##  [1] -4.5 -3.5 -2.5 -1.5 -0.5  0.5  1.5  2.5  3.5  4.5
## 
## $xname
## [1] "g"
## 
## $equidist
## [1] TRUE
## 
## attr(,"class")
## [1] "histogram"
g <- rtns[,3]
h<-hist(g, breaks=10, density=10, col="lightgray", xlab="Accuracy", main="Overall") 
    xfit<-seq(min(g),max(g),length=40) 
    yfit<-dnorm(xfit,mean=mean(g),sd=sd(g)) 
    yfit <- yfit*diff(h$mids[1:2])*length(g) 
    lines(xfit, yfit, col="black", lwd=2)

h
## $breaks
##  [1] -6 -5 -4 -3 -2 -1  0  1  2  3  4  5  6
## 
## $counts
##  [1]    1    3   13  132  423 1072 1259  464   90    7    1    1
## 
## $density
##  [1] 0.0002885170 0.0008655511 0.0037507213 0.0380842470 0.1220427005
##  [6] 0.3092902481 0.3632429313 0.1338718984 0.0259665320 0.0020196192
## [11] 0.0002885170 0.0002885170
## 
## $mids
##  [1] -5.5 -4.5 -3.5 -2.5 -1.5 -0.5  0.5  1.5  2.5  3.5  4.5  5.5
## 
## $xname
## [1] "g"
## 
## $equidist
## [1] TRUE
## 
## attr(,"class")
## [1] "histogram"
g <- rtns[,4]
h<-hist(g, breaks=10, density=10, col="lightgray", xlab="Accuracy", main="Overall") 
    xfit<-seq(min(g),max(g),length=40) 
    yfit<-dnorm(xfit,mean=mean(g),sd=sd(g)) 
    yfit <- yfit*diff(h$mids[1:2])*length(g) 
    lines(xfit, yfit, col="black", lwd=2)

h
## $breaks
##  [1] -5 -4 -3 -2 -1  0  1  2  3  4  5
## 
## $counts
##  [1]    3   19  119  431 1051 1258  478   88   16    3
## 
## $density
##  [1] 0.0008655511 0.0054818234 0.0343335257 0.1243508367 0.3032313907
##  [6] 0.3629544143 0.1379111368 0.0253894980 0.0046162724 0.0008655511
## 
## $mids
##  [1] -4.5 -3.5 -2.5 -1.5 -0.5  0.5  1.5  2.5  3.5  4.5
## 
## $xname
## [1] "g"
## 
## $equidist
## [1] TRUE
## 
## attr(,"class")
## [1] "histogram"
g <- rtns[,5]
h<-hist(g, breaks=10, density=10, col="lightgray", xlab="Accuracy", main="Overall") 
    xfit<-seq(min(g),max(g),length=40) 
    yfit<-dnorm(xfit,mean=mean(g),sd=sd(g)) 
    yfit <- yfit*diff(h$mids[1:2])*length(g) 
    lines(xfit, yfit, col="black", lwd=2)

h
## $breaks
##  [1] -5 -4 -3 -2 -1  0  1  2  3  4  5
## 
## $counts
##  [1]    2   20  122  422 1118 1205  479   82   10    6
## 
## $density
##  [1] 0.000577034 0.005770340 0.035199077 0.121754183 0.322562031
##  [6] 0.347663012 0.138199654 0.023658396 0.002885170 0.001731102
## 
## $mids
##  [1] -4.5 -3.5 -2.5 -1.5 -0.5  0.5  1.5  2.5  3.5  4.5
## 
## $xname
## [1] "g"
## 
## $equidist
## [1] TRUE
## 
## attr(,"class")
## [1] "histogram"
g <- rtns[,6]
h<-hist(g, breaks=10, density=10, col="lightgray", xlab="Accuracy", main="Overall") 
    xfit<-seq(min(g),max(g),length=40) 
    yfit<-dnorm(xfit,mean=mean(g),sd=sd(g)) 
    yfit <- yfit*diff(h$mids[1:2])*length(g) 
    lines(xfit, yfit, col="black", lwd=2)

h
## $breaks
##  [1] -7 -6 -5 -4 -3 -2 -1  0  1  2  3  4  5  6
## 
## $counts
##  [1]    1    1   10   42  112  424 1099 1166  460  123   22    4    2
## 
## $density
##  [1] 0.000288517 0.000288517 0.002885170 0.012117715 0.032313907
##  [6] 0.122331218 0.317080208 0.336410848 0.132717830 0.035487594
## [11] 0.006347374 0.001154068 0.000577034
## 
## $mids
##  [1] -6.5 -5.5 -4.5 -3.5 -2.5 -1.5 -0.5  0.5  1.5  2.5  3.5  4.5  5.5
## 
## $xname
## [1] "g"
## 
## $equidist
## [1] TRUE
## 
## attr(,"class")
## [1] "histogram"
g <- rtns[,7]
h<-hist(g, breaks=10, density=10, col="lightgray", xlab="Accuracy", main="Overall") 
    xfit<-seq(min(g),max(g),length=40) 
    yfit<-dnorm(xfit,mean=mean(g),sd=sd(g)) 
    yfit <- yfit*diff(h$mids[1:2])*length(g) 
    lines(xfit, yfit, col="black", lwd=2)

h
## $breaks
##  [1] -5 -4 -3 -2 -1  0  1  2  3  4  5
## 
## $counts
##  [1]    3    3   65  373 1228 1346  396   48    3    1
## 
## $density
##  [1] 0.0008655511 0.0008655511 0.0187536065 0.1076168494 0.3542989036
##  [6] 0.3883439123 0.1142527409 0.0138488171 0.0008655511 0.0002885170
## 
## $mids
##  [1] -4.5 -3.5 -2.5 -1.5 -0.5  0.5  1.5  2.5  3.5  4.5
## 
## $xname
## [1] "g"
## 
## $equidist
## [1] TRUE
## 
## attr(,"class")
## [1] "histogram"
m <- pobs(as.matrix(cbind(rtns)))

splom2(m)

Normal copula

n.cop <- normalCopula(dim=7, dispstr = 'un')
set.seed(500)
fit <- fitCopula(n.cop,m,method='ml')
fit
## fitCopula() estimation based on 'maximum likelihood'
## and a sample of size 3466.
##          Estimate Std. Error  z value Pr(>|z|)    
## rho.1   9.714e-01  6.833e-04 1421.637  < 2e-16 ***
## rho.2   4.661e-01  1.195e-02   38.994  < 2e-16 ***
## rho.3  -3.060e-02  1.702e-02   -1.798  0.07225 .  
## rho.4  -3.929e-02  1.701e-02   -2.310  0.02092 *  
## rho.5  -1.696e-02  1.705e-02   -0.995  0.31989    
## rho.6  -1.430e-02  1.705e-02   -0.839  0.40166    
## rho.7   4.915e-01  1.157e-02   42.464  < 2e-16 ***
## rho.8  -2.463e-02  1.703e-02   -1.447  0.14801    
## rho.9  -4.021e-02  1.701e-02   -2.364  0.01809 *  
## rho.10 -1.877e-02  1.705e-02   -1.101  0.27075    
## rho.11 -1.387e-02  1.705e-02   -0.814  0.41593    
## rho.12 -4.980e-02  1.699e-02   -2.931  0.00338 ** 
## rho.13  1.646e-06  1.705e-02    0.000  0.99992    
## rho.14  1.580e-03  1.705e-02    0.093  0.92620    
## rho.15 -5.006e-03  1.706e-02   -0.293  0.76914    
## rho.16 -1.738e-02  1.705e-02   -1.019  0.30805    
## rho.17  1.913e-02  1.705e-02    1.122  0.26171    
## rho.18  8.482e-03  1.705e-02    0.497  0.61893    
## rho.19  9.321e-03  1.705e-02    0.547  0.58467    
## rho.20 -1.856e-02  1.705e-02   -1.089  0.27628    
## rho.21  1.395e-02  1.705e-02    0.818  0.41315    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## The maximized loglikelihood is  5474 
## Optimization converged
## Number of loglikelihood evaluations:
## function gradient 
##      196       28
rho <-coef(fit)
copula_dist <- mvdc(copula=normalCopula(rho,dim=7,dispstr = 'un'), 
                    margins=c("norm","norm","norm","norm","norm","norm","norm"),
                    paramMargins=list(list(mean=mean(rtns[,1]), sd=sd(rtns[,1])),
                                      list(mean=mean(rtns[,2]), sd=sd(rtns[,2])),
                                      list(mean=mean(rtns[,3]), sd=sd(rtns[,3])),
                                      list(mean=mean(rtns[,4]), sd=sd(rtns[,4])),
                                      list(mean=mean(rtns[,5]), sd=sd(rtns[,5])),
                                      list(mean=mean(rtns[,6]), sd=sd(rtns[,6])),
                                      list(mean=mean(rtns[,7]), sd=sd(rtns[,7]))
                                      )
                    )

sim <- rmvdc(copula_dist, 3466)
## Warning: 'rmvdc' is deprecated.
## Use 'rMvdc' instead.
## See help("Deprecated")
plot(rtns[,1],rtns[,2],main='normal-copula DJIA-GSPC')
points(sim[,1],sim[,2],col='purple')
legend('bottomright',c('Observed','Simulated'),col=c('black','purple'),pch=21)

plot(rtns[,1],rtns[,3],main='normal-copula DJIA-NDX')
points(sim[,1],sim[,3],col='purple')
legend('bottomright',c('Observed','Simulated'),col=c('black','purple'),pch=21)

plot(rtns[,1],rtns[,4],main='normal-copula DJIA-GDAXI')
points(sim[,1],sim[,4],col='purple')
legend('bottomright',c('Observed','Simulated'),col=c('black','purple'),pch=21)

plot(rtns[,1],rtns[,5],main='normal-copula DJIA-FCHI')
points(sim[,1],sim[,5],col='purple')
legend('bottomright',c('Observed','Simulated'),col=c('black','purple'),pch=21)

plot(rtns[,1],rtns[,6],main='normal-copula DJIA-SSEC')
points(sim[,1],sim[,6],col='purple')
legend('bottomright',c('Observed','Simulated'),col=c('black','purple'),pch=21)

plot(rtns[,1],rtns[,7],main='normal-copula DJIA-SENSEX')
points(sim[,1],sim[,7],col='purple')
legend('bottomright',c('Observed','Simulated'),col=c('black','purple'),pch=21)

t distr copula

t.cop <- tCopula(dim=7, dispstr = 'un')
set.seed(500)
fit <- fitCopula(t.cop,m,method='ml')
fit
## fitCopula() estimation based on 'maximum likelihood'
## and a sample of size 3466.
##          Estimate Std. Error  z value Pr(>|z|)    
## rho.1   0.9713667  0.0007552 1286.156  < 2e-16 ***
## rho.2   0.5145676  0.0124411   41.360  < 2e-16 ***
## rho.3  -0.0293183  0.0177551   -1.651  0.09869 .  
## rho.4  -0.0316248  0.0178154   -1.775  0.07588 .  
## rho.5  -0.0133095  0.0177248   -0.751  0.45272    
## rho.6  -0.0034288  0.0176745   -0.194  0.84618    
## rho.7   0.5414400  0.0120706   44.856  < 2e-16 ***
## rho.8  -0.0242672  0.0177692   -1.366  0.17204    
## rho.9  -0.0333458  0.0178330   -1.870  0.06150 .  
## rho.10 -0.0159670  0.0177309   -0.901  0.36785    
## rho.11 -0.0038308  0.0176985   -0.216  0.82864    
## rho.12 -0.0473888  0.0177987   -2.662  0.00776 ** 
## rho.13 -0.0045290  0.0179107   -0.253  0.80037    
## rho.14 -0.0062200  0.0177792   -0.350  0.72645    
## rho.15  0.0038222  0.0177604    0.215  0.82960    
## rho.16 -0.0198613  0.0175692   -1.130  0.25828    
## rho.17  0.0191351  0.0175063    1.093  0.27438    
## rho.18  0.0013937  0.0175002    0.080  0.93653    
## rho.19  0.0018595  0.0175989    0.106  0.91585    
## rho.20 -0.0199084  0.0174932   -1.138  0.25509    
## rho.21  0.0196046  0.0174675    1.122  0.26172    
## df     10.5921609  0.6350647   16.679  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## The maximized loglikelihood is  5681 
## Optimization converged
## Number of loglikelihood evaluations:
## function gradient 
##      184       33
rho <- coef(fit)[1:21]
df <- coef(fit)[22]
copula_dist <- mvdc(copula=normalCopula(rho,dim=7,dispstr = 'un'), 
                    margins=c("t","t","t","t","t","t","t"),
                    paramMargins=list(list(df = df),
                                      list(df = df),
                                      list(df = df),
                                      list(df = df),
                                      list(df = df),
                                      list(df = df),
                                      list(df = df)
                                      )
                    )

sim <- rmvdc(copula_dist, 3466)
## Warning: 'rmvdc' is deprecated.
## Use 'rMvdc' instead.
## See help("Deprecated")
plot(rtns[,1],rtns[,2],main='t-copula DJIA-GSPC')
points(sim[,1],sim[,2],col='purple')
legend('bottomright',c('Observed','Simulated'),col=c('black','purple'),pch=21)

plot(rtns[,1],rtns[,3],main='t-copula DJIA-NDX')
points(sim[,1],sim[,3],col='purple')
legend('bottomright',c('Observed','Simulated'),col=c('black','purple'),pch=21)

plot(rtns[,1],rtns[,4],main='t-copula DJIA-GDAXI')
points(sim[,1],sim[,4],col='purple')
legend('bottomright',c('Observed','Simulated'),col=c('black','purple'),pch=21)

plot(rtns[,1],rtns[,5],main='t-copula DJIA-FCHI')
points(sim[,1],sim[,5],col='purple')
legend('bottomright',c('Observed','Simulated'),col=c('black','purple'),pch=21)

plot(rtns[,1],rtns[,6],main='t-copula DJIA-SSEC')
points(sim[,1],sim[,6],col='purple')
legend('bottomright',c('Observed','Simulated'),col=c('black','purple'),pch=21)

plot(rtns[,1],rtns[,7],main='t-copula DJIA-SENSEX')
points(sim[,1],sim[,7],col='purple')
legend('bottomright',c('Observed','Simulated'),col=c('black','purple'),pch=21)

Clayton copula

c.cop <- claytonCopula(dim=7)
set.seed(500)

fit <- fitCopula(c.cop,m,method='ml')
fit
## fitCopula() estimation based on 'maximum likelihood'
## and a sample of size 3466.
##       Estimate Std. Error z value Pr(>|z|)    
## param 0.095633   0.005807   16.47   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## The maximized loglikelihood is  226.4 
## Optimization converged
## Number of loglikelihood evaluations:
## function gradient 
##       24        6
theta <- coef(fit)[1]
copula_dist <- mvdc(copula=claytonCopula(theta,dim=7), 
                    margins=c("norm","norm","norm","norm","norm","norm","norm"),
                    paramMargins=list(list(mean=mean(rtns[,1]), sd=sd(rtns[,1])),
                                      list(mean=mean(rtns[,2]), sd=sd(rtns[,2])),
                                      list(mean=mean(rtns[,3]), sd=sd(rtns[,3])),
                                      list(mean=mean(rtns[,4]), sd=sd(rtns[,4])),
                                      list(mean=mean(rtns[,5]), sd=sd(rtns[,5])),
                                      list(mean=mean(rtns[,6]), sd=sd(rtns[,6])),
                                      list(mean=mean(rtns[,7]), sd=sd(rtns[,7]))
                                      )
                    )

sim <- rmvdc(copula_dist, 3466)
## Warning: 'rmvdc' is deprecated.
## Use 'rMvdc' instead.
## See help("Deprecated")
plot(rtns[,1],rtns[,2],main='Clayton-copula DJIA-GSPC')
points(sim[,1],sim[,2],col='purple')
legend('bottomright',c('Observed','Simulated'),col=c('black','purple'),pch=21)

plot(rtns[,1],rtns[,3],main='Clayton-copula DJIA-NDX')
points(sim[,1],sim[,3],col='purple')
legend('bottomright',c('Observed','Simulated'),col=c('black','purple'),pch=21)

plot(rtns[,1],rtns[,4],main='Clayton-copula DJIA-GDAXI')
points(sim[,1],sim[,4],col='purple')
legend('bottomright',c('Observed','Simulated'),col=c('black','purple'),pch=21)

plot(rtns[,1],rtns[,5],main='Clayton-copula DJIA-FCHI')
points(sim[,1],sim[,5],col='purple')
legend('bottomright',c('Observed','Simulated'),col=c('black','purple'),pch=21)

plot(rtns[,1],rtns[,6],main='Clayton-copula DJIA-SSEC')
points(sim[,1],sim[,6],col='purple')
legend('bottomright',c('Observed','Simulated'),col=c('black','purple'),pch=21)

plot(rtns[,1],rtns[,7],main='Clayton-copula DJIA-SENSEX')
points(sim[,1],sim[,7],col='purple')
legend('bottomright',c('Observed','Simulated'),col=c('black','purple'),pch=21)

Gumbel copula

c.cop <- gumbelCopula(dim=7)
set.seed(500)

fit <- fitCopula(c.cop,m,method='ml')
## Warning in .local(copula, tau, ...): tau is out of the range [0, 1]

## Warning in .local(copula, tau, ...): tau is out of the range [0, 1]

## Warning in .local(copula, tau, ...): tau is out of the range [0, 1]

## Warning in .local(copula, tau, ...): tau is out of the range [0, 1]

## Warning in .local(copula, tau, ...): tau is out of the range [0, 1]

## Warning in .local(copula, tau, ...): tau is out of the range [0, 1]

## Warning in .local(copula, tau, ...): tau is out of the range [0, 1]

## Warning in .local(copula, tau, ...): tau is out of the range [0, 1]

## Warning in .local(copula, tau, ...): tau is out of the range [0, 1]

## Warning in .local(copula, tau, ...): tau is out of the range [0, 1]

## Warning in .local(copula, tau, ...): tau is out of the range [0, 1]

## Warning in .local(copula, tau, ...): tau is out of the range [0, 1]

## Warning in .local(copula, tau, ...): tau is out of the range [0, 1]

## Warning in .local(copula, tau, ...): tau is out of the range [0, 1]

## Warning in .local(copula, tau, ...): tau is out of the range [0, 1]

## Warning in .local(copula, tau, ...): tau is out of the range [0, 1]
fit
## fitCopula() estimation based on 'maximum likelihood'
## and a sample of size 3466.
##       Estimate Std. Error z value Pr(>|z|)    
## param 1.043675   0.004049   257.7   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## The maximized loglikelihood is  89.32 
## Optimization converged
## Number of loglikelihood evaluations:
## function gradient 
##       39       10
theta <- coef(fit)[1]
copula_dist <- mvdc(copula=gumbelCopula(theta,dim=7), 
                    margins=c("norm","norm","norm","norm","norm","norm","norm"),
                    paramMargins=list(list(mean=mean(rtns[,1]), sd=sd(rtns[,1])),
                                      list(mean=mean(rtns[,2]), sd=sd(rtns[,2])),
                                      list(mean=mean(rtns[,3]), sd=sd(rtns[,3])),
                                      list(mean=mean(rtns[,4]), sd=sd(rtns[,4])),
                                      list(mean=mean(rtns[,5]), sd=sd(rtns[,5])),
                                      list(mean=mean(rtns[,6]), sd=sd(rtns[,6])),
                                      list(mean=mean(rtns[,7]), sd=sd(rtns[,7]))
                                      )
                    )

sim <- rmvdc(copula_dist, 3466)
## Warning: 'rmvdc' is deprecated.
## Use 'rMvdc' instead.
## See help("Deprecated")
plot(rtns[,1],rtns[,2],main='Gumbel-copula DJIA-GSPC')
points(sim[,1],sim[,2],col='purple')
legend('bottomright',c('Observed','Simulated'),col=c('black','purple'),pch=21)

plot(rtns[,1],rtns[,3],main='Gumbel-copula DJIA-NDX')
points(sim[,1],sim[,3],col='purple')
legend('bottomright',c('Observed','Simulated'),col=c('black','purple'),pch=21)

plot(rtns[,1],rtns[,4],main='Gumbel-copula DJIA-GDAXI')
points(sim[,1],sim[,4],col='purple')
legend('bottomright',c('Observed','Simulated'),col=c('black','purple'),pch=21)

plot(rtns[,1],rtns[,5],main='Gumbel-copula DJIA-FCHI')
points(sim[,1],sim[,5],col='purple')
legend('bottomright',c('Observed','Simulated'),col=c('black','purple'),pch=21)

plot(rtns[,1],rtns[,6],main='Gumbel-copula DJIA-SSEC')
points(sim[,1],sim[,6],col='purple')
legend('bottomright',c('Observed','Simulated'),col=c('black','purple'),pch=21)

plot(rtns[,1],rtns[,7],main='Gumbel-copula DJIA-SENSEX')
points(sim[,1],sim[,7],col='purple')
legend('bottomright',c('Observed','Simulated'),col=c('black','purple'),pch=21)